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Abstract 

The quark number susceptibihty is expressed as an integral in terms of dressed quark propagator 
and dressed vector vertex. It is then investigated with the latter two- and three-point functions 
confronted with a Dyson-Schwinger equation model which accommodates both finite temperature 
and baryon chemical potential. The critical end point in the phase diagram is identified and the 
behavior of the quark number susceptibility around the critical end point is highlighted. The 
qualitative features found agree with recent lattice QCD simulation results. 
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I. INTRODUCTION 



Underlying our understanding of strong interaction, Quantum Chromodynamics (QCD) 
is remarkable in that its fundamental degrees of freedom - quarks and gluons are confined 
and the (approximate) chiral symmetry in the light quark sector is dynamically broken. On 
the other hand, at sufficiently high temperatures and/or baryon chemical potentials, QCD 
is believed to undergo a phase transition into the deconfined quark gluon plasma (QGP) 
(for a recent review, see, e.g., [iQ, and probably at the same time, the spontaneously broken 
chiral symmetry gets restored ^. It is the main task of the on-going BNL-RHIC and one 
of the goals of CERN-LHC to create and probe such a new state of QCD matter. 

The phase transition of confined hadronic matter to QGP at finite temperature and 
vanishing baryon chemical potential has been extensively studied by lattice QCD simulations 



as well as phenomenological models. The chiral 
in the chiral limit is likely to be of second order [3 



case of physical finite current quark masses 



m 



lase transition of QCD with two fiavors 
and turns to a smooth crossover in the 
6|. At finite chemical potentials and lower 



temperatures, more and more robust lattice simulations indicate that QCD experiences a 



first-order phase transition 0, H]. Such a scenario of QCD phase transition suggests that 
there should exist a critical end point (CEP) on the phase diagram - the point where the 
first-order phase transition line terminates (for a recent review on CEP, and especially the 
argument that CEP must exist, see [9l). 

The search for the QCD CEP has attracted considerable attention over the years from 
both theoretical and experimental aspects Q, Q, Q, Q, Q, Q, Q- Note that the end 
point of a first-order phase transition line is a critical point of the second order - in this 
connection, the water liquid-vapor phase transition makes a more familiar example. This 
means that the CEP would be characterized by enhanced long wave-length fiuctuations 
which lead to singularity in thermodynamic observables. The authors of Ref. |10j] suggest 
that quark number susceptibility should develop some singularity near the CEP and its 



divergence could be a sig nal of the CEP. From then on, a lot of p 



investigations 
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nenomeno 



22| as well as lattice simulations 
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ogical model 



25| toward this 



susceptibility have been attempted, all aimed at determining the location of the CEP on the 
temperature (T)-chemical potential (/i) plane. On the experimental side, recent plans for 
verification of the CEP in heavy ion collisions have concentrated on energy scans to access 
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as broad as possible range of T and /i values necessary for the search for the CEP. It is 
argued that the range of T and fi values for the CEP search can be accessible by combined 



results from energy scans at FAIR, SPS and RHIC [16|]. Also notable is that, apart from 
susceptibilities, transport coefficients, like the ratio of viscosity to entropy density {i]/s) can 
be other probes for the CEP 1^, 16, 26|. 

In the present work, we extend our previous continuum study of quark number suscep- 
tibility to finite chemical potential in a Dyson- Schwinger equation model and aim at 
locating the CEP. We first derive a closed integral expression for the susceptibility in terms 
of QCD dressed quark propagator and dressed vector vertex. This integral formula proves 
to reproduce the exact analytical result in the case of a free quark gas. In the next sec- 
tion, with a model gluon propagator input, the 2- and 3-point functions are confronted with 
Dyson- Schwinger equation (DSE) and Bethe-Salpeter equation (BSE). With these results, 
the quark number susceptibility is then calculated and the characteristic behavior of the 
susceptibility near the CEP is highlighted and its implications analyzed. The last section is 
devoted to some concluding remarks. Throughout this work, we work with Euclidean space 
metric: {7^,71.} = 26f,„. 

II. ANALYTICAL RESULTS 

A. Derivation of a closed integral expression for the quark number susceptibility 

The quark number susceptibility is defined as the derivative of quark number density 
with respect to the baryon chemical potential 

-I- « 

To start our derivation, we write the quark number density as 

j^tr^[G{p, (2) 

where G denotes the quark propagator, Nc, Nf the color factor and the number of fiavors 
respectively, and the trace is taken over Dirac indices. Here, the minus sign comes from the 
anti-commutation property of the fermion fields. 

Substituting Eq. (2) into Eq. (1) and adopting the identity 

dG{p,fi) dG-\p,fi) 

— — = (-)G'(p,^) ^ G{p,fi), (3) 



we have 



Recall that the well-known Ward identity reads 

where p denotes the relative momentum of the vector vertex and the corresponding total 
momentum vanishes. Note that at finite temperature and chemical potential, the fourth 
component of momentum p4 = cun + il-i, with cun being fermion Matsubara frequencies. Then 

= (6) 

Putting this equation into Eq. (4) and replacing the integration over the fourth component 
of momentum with explicit summation over Matsubara frequencies, we get the quark number 
susceptibility at finite temperature and chemical potential in the imaginary-time formalism 

xiT,f,)^{-)N,NfT J2 I 7A-,tr,[G{pn)r,{p^,0)G{pM, (7) 

n=— oo ^ ' 

where p„ = ( p*', -|- ifj) with fermion frequencies a;„ = {2n + 1)ttT. 

Thus a model-independent closed integral formula is obtained, which expresses the quark 
number susceptibility in terms of the dressed quark propagator and the dressed vector vertex, 
both of the latter objects being basic quantities in quantum field theory. The DSE-BSE 
approach provides an ideal framework to confront these quantities non-perturbatively and 
then calculate the quark number susceptibility. 



B. The free quark number susceptibility 

Before going into the model calculation of the quark number susceptibility, let us apply 
the foregoing analytical expression for the susceptibility to the case of a free quark gas. In 
this case, the dressed vector vertex reduces to the bare one 

r4(Pn, 0)^74 (8) 

and the free quark propagator reads 



with m the degenerated hght quark current mass. Substituting Eqs. (8) and (9) into Eq. 
(7) and evaluating the trace, one gets 

where E^; = \/f^ + m?. With the aid of the identity 

y — - — = ^ X (11) 



2C(J„ — X 

n=— oo 



one immediately arrives at 



This equation can be rewritten as 

X 



This integral is nothing but the free quark gas density from the Fermi-Dirac statistics. Eq. 
(13) is just a re-expression of the definition of the quark number susceptibility (Eq.(l)) for 
a free quark gas. The integration of Eq. (13) can be analytically solved in the chiral limit 
(m = 0) 

X^^^^(T,/i)=iV^(T2 + ^). (14) 

We see that the integral formula Eq. (7) for the quark number susceptibility reproduces 
exactly the Fermi-Dirac statistics result for a free quark gas. This is to be expected in 
advance and can be regarded as a consistent check of Eq. (7). 



III. MODEL CALCULATION 

A. The dressed quark propagator 

In order to calculate the quark number susceptibility, we have to evaluate first the dressed 
quark propagator and dressed vector vertex at finite temperature and chemical potential. 
This will be conducted systematically and consistently in the rainbow-DSE and ladder-BSE 
framework with a separable effective interaction model. 

The DSE-BSE provides a non-perturbative continuum approach for the exploration of 
strong interaction physics. Its practical form, namely the rainbow-ladder truncation, has 

5 



given us fascinating insights into the remarkable dynamical chiral symmetry breaking and 



confinement 



28 



291] and found successful applications in describing light mesons in the 



pseudo-scalar and vector channels (for review articles, see 30|]). The rainbow DSE for the 
quark propagator reads 



G{p) 



-1 



i'y ■ p + m + - 



(2^ 



g^D%f{p-q)^,G{q)^, 



and the ladder inhomogeneous BSE for the dressed vector vertex reads 

4 r d^q 



9'D;fJ{p - q)%G{q^)r^iq, P)G{q^h^ 



(15) 



(16) 



3 J (27r)4 

where P denotes the total momentum, q± = q ± P/2 and D^JJ{p — q) the effective gluon 
propagator, which is usually a phenomenological input in practice (There are also some 
attempts to explore the analytical structure of the gluon propagator from numerical solutions 



of coupled DSEs of quarks, gluons and ghosts and from fits to lattice date, see Ref. 
and references therein") . In Refs. 



31 



29| 



32] , the authors proposed a confining, separable model 



gluon propagator with the rank-1 form 



g'DlfJ{p~q) = 6,Mo{p')fo{q' 



(17) 



where /o(p^) = exp(— p^/A^). This model is found to be successful in describing light fiavor 
pseudo-scalar and vector meson observables with parameters Aq = 0.678 GeV, DqA^ = 128.0 
and the degenerated light quark mass m = 6.6 MeV 32]. 

In recent years, the rainbow-ladder DSE models were extended to finite teniperature 
and chemical potential to investigate QCD thermodynamics (for a review, see 33|]). As 
for the rank-1 confining separable model, the extension to finite temperature and chemical 
potential is accomplished by transcription of the Euclidean quark four-momentum via q — > 



In = (^, i^n), with LJn = + ifi, and uo ucw parameters are introduced [32|, 13^, l35|. This 



means the effective gluon propagator at finite temperature and chemical potential is modeled 

as 



g'D';/J{pk - qn) = S.MoipDfoiq'n)- 



(18) 



On the other hand, the finite temperature and chemical potential version of the rainbow-DSE 
of quark propagator (Eq. 15) reads 



G \pk) = Pk + m + -T I 



cP'q 
(2^ 



(19) 



and the quark propagator is generally decomposed as 

G-\pk) = ^l ■ PA{PI) + ^l4^kC{pl) + B{pI). (20) 

The rank-1 separable model leads to the rainbow-DSE solution: A{pl) = C{pl) = 1 and 
B{pf.) = m + b(T, fi)fo(ji^), where b{T,fi) satisfies the following equation 

b{T,,) - ^D.Tl^J + + (21) 

Of course, that the vector self-energy amplitudes A and C are unity does not cover all the 
physical effects, but this respects the nature of the separable model - the correspondence 
between quark pro pagator scalar self-energy amplitude B and the separable gluon propagator 
form factor /o [35|]- Furthermore, it is an open question how the chemical potential /i enters 
the gluon propagator and the forgoing transcription of the fourth-component momentum 
P4 — > ujn = ojn + i/i has uot bccu justified. However, this is again the requirement of the 
nature of the separable model and we think it worthwhile to allow such a uncomplicated 
extension to finite T and /i in which the Matsubara modes are not coupled and one could sum 
them as easily as possible and then make a first survey of the thermodynamic observables 



35l |. Fortunately the results turn out to be satisfac tory - the correspondence nature of the 



lail 



separable model plays a vital and correct role here 

Eq. (21) is numerically solved. Note that this simple model quite facilitates summa- 
tion over fermion Matsubara frequencies as well as the 3-momentum integration due to its 
ultraviolet finiteness characterized by the Gaussian function fo{q^). As a result, no renor- 
malization is needed. The imaginary part of b{T, /i) is found to be very small (of order 
~ 10"^*^) for all temperatures and chemical potentials while the real part takes strikingly 
different values at different T and fi regions, hence serving as the order parameter. Fig. 1 
shows the /^-dependent behavior of Re6(T, /i) at several fixed temperatures. One sees that at 
low temperatures, Reb(T, jj) demonstrates an abrupt discontinuous decrease to a very small 
value at some critical chemical potential, which is typical of a phase transition of first order 



and resembles the observations of Nambu-Jona-Lasinio (NJL) model 36|]. Nevertheless, at 
higher temperatures, Re6(T, /i) shows a smooth, gradual decrease as /x increases, implying 
an analytic crossover. So one can naturally expect a CEP at which the first-order phase 
transition and the crossover connect. This will be verified by determining the most singular 
point of the quark number susceptibility, which is illustrated in subsection C. To summarize 



this subsection, we point out that the observations of Fig. 1 agree with the scenario of QCD 
phase transition shortly reviewed in the introduction. 
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FIG. 1: //-dependence of Iieb{T,fi) at three fixed temperatures. The evolution from a first-order 
phase transition to a crossover is shown. The different monotonous behavior of Re6(T, fj,) in the 
symmetry broken phase for T = 90 MeV and T = 110 MeV may be ascribed to the artifact of the 
separable model. 



B. The dressed vector vertex 

We now turn to the dressed vector vertex. The finite temperature and chemical po- 
tential version of the inhomogeneous ladder BSE for the dressed vector vertex (the fourth 
component) Eq. (16) with vanishing total momentum reads 

Uipk, 0) = 74 - 3T 5^ / —^g^DfJipk - qn)lpG{qn)V,{qn, 0)G(g„)7.. (22) 

n=— 00 ^ ' 

where the model gluon propagator g'^Dpl^{pk — qn) is given by Eq. (18). 

From Lorentz structure analysis, the dressed vector vertex with vanishing total momen- 



tum at T = = /i can be decomposed as [SOj 



r^(p, 0) = ai(p^)7/. + "2(^)7 ■ PPm - "3(/)^Pm- (23) 



This can also be easily seen from the Ward identity iT{p, 0) = dG ^{p)/dp^ with G ^{p) = 
27 ■ pAlp"^) + B{p^). In Eq. (23), ai{p^) is the dominant term, while the other two terms 



contribute so little that they can even be neglected 30(]. At finite temperature and chemical 



potential, we keep the corresponding "leading" term in V^{pi^,{]) and add another term 

r4(Pfc, 0) = a{pi)^, - iPipi)uk. (24) 

Here we stress that the second term is vital and indispensable for describing the QCD phase 
transition at finite temperature and chemical potential. Some explanation goes as follows 
in this connection. From the Ward identity Eq. (6) one sees that this term arises from the 
derivative of the scalar self-energy amplitude -B(p|) of the inverse dressed quark propagator 
Eq. (20) with respect to the chemical potential fi. Note that -B(pf,) experiences an abrupt 
decrease at some critical fi when a true phase transition occurs (on the contrary, the vector 



self-energy amplitudes A(j)l) and C(p|) evolve smoothly with T and /i [32]) and hence makes 
a considerable contribution to T4^{pk,0) at that point. 

Putting the ansatz for the dressed vector vertex Eq. (24), the dressed quark propagator 
Eq. (20) and the confining separable model gluon propagator Eq. (18) into Eq. (22), one 
can transform the ladder BSE into two coupled non-linear equations 



" ~^ " ,hJ (^^°** ' |,7 + (m + (,/„(5?))T ^ 

""^ J (2^)3 ) [5f+(m + Vo(5?))T ^ 



with aipl) = 1 — u(T, fi)fQ{pl) and I3{p\)ujk = v(T, fi)fQ(j)l). With b(T,fi) obtained from 
Eq. (21), these two equations are numerically solved. Here, one also notes that the Gaussian 
function /o($f) appearing in the numerators of the integrands justifies a numerical cut-off 
for summation over Matsubara frequencies as well as 3-D momentum integration. 



C. Results for the quark number susceptibility 

Having completed the calculation of the dressed quark propagator and the dressed vector 
vertex, we are now in a position to calculate the quark number susceptibility. The numerical 
results are shown in Fig. 2, where x{T, /i) is normalized by the free quark quark number 
susceptibility in the chiral limit x^^'^^{T, 11) (Eq. (14)) and hence dimensionless. 
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FIG. 2: T-dependence of x(^) ^) &t four fixed chemical potentials. From left to right, /i = 
179, 164, 149, 134 MeV, respectively. The CEP is identified as the diverging cusp developed in 
the solid line at Tcep = 117 MeV,/icEP = 164 MeV. 

Let us analyse the observations from Fig. 2. The most singular point develops along 
the /icEP = 164 MeV solid line at Tcep = 117 MeV. This is the CEP, at which the quark 
number susceptibility manifests itself as a diverging cusp that should translate into event 



by event fluctuations of baryon multiplicity 



22l |. For larger /i = 179 MeV (the dashed line), 



x{T, fi) also jumps across the chiral phase transition and has a discontinuity due to its first- 
order character. The peak is still pronounced but its height much suppressed. For smaller 
chemical potentials (/i = 149, 134 MeV, the dash-dotted and the dotted line), one enters a 
crossover region, in which the discontinuity of x{T, f^) at the transition vanishes and x{T, f^) 
evolves gradually with T in a smoother and smoother manner. 

Another observation from Fig. 2 is that far below the chiral phase transition, x(T, /i) 
is greatly suppressed and at asymptotic temperatures, x(T, /i) tends unambiguously to the 
value of susceptibility of a free quark gas. In the latter respect, our demonstration is more 
akin to that of lattice predictions [23|, l25|] than NJL-type models [l8|, l20| , where x decreases 
too fast at asymptotic temperatures or chemical potentials. This may be accounted for by the 
fact that the momentum-dependent interaction DSE model represents a better asymptotic 
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freedom behavior. 

Our calculation determines the CEP location at Tqep = 117 MeV and /icEP = 164 MeV. 
The most recent experimental estimate by extracting 77/5 as a function of T and fi from an 
elliptic flow excitation function gives Tqep ~ 165 — 170 MeV and /icEP ~ 150 — 180 MeV 



Iq . which is in good agreement with lattice prediction with almost realistic pion mass 24 1. 



So very encouragingly, our results for /icEP lies exactly in this interval, whereas our Tcep 
deviates a lot, for the pseudo-critical temperature Tc in our calculation takes a relatively 
lower value (see below). As for NJL-type models, they predict a much higher /xcep isl- I20I. 



221]. In fact, NJL-type models have a common tendency to lead to the CEP at relatively 



high /icEP ~ 300 MeV [22 1. 

Having identified the CEP, we now show the complete phase diagram obtained from our 
model, which assumes the shape sketched in Fig. 3. Two separate phases are shown in 



160 I I I I I I I I I I I I I I I — ■—I—' — r-"— 1 — i I i I — '—I—' — |— ^ 
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FIG. 3: The phase diagram characterized by Re6(T, ^). The solid line represents the first-order 
phase transition and the dotted line the crossover transition. The filled circle indicates the CEP. 

the T — fi plane: the chiral symmetry restored phase and the phase with (approximate) 
chiral symmetry dynamically broken. The solid line represents a first-order phase transition 
line, which begins from fic = 286 MeV at vanishing temperature and ends at the CEP 
(Tcep = 117 MeV, /xcep = 164 MeV). Beyond the Tcep, one ends up with a crossover 
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(the dotted line) with a pseudo-critical temperature Tc — 150 MeV at vanishing chemical 
potential. 

Last but not least, to check the parameter sensitivity of the location of the critical end 
point, we first vary the degenerated current quark mass m, since the value of m is believed 
to be closely associated with the nature of the phase transition. The numerical results is 
shown in Fig. 4. 
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FIG. 4: The location of the critical end point for m = 5.2 MeV and ra = 8.0 MeV with "physical" 
Dq = 128.0/(0.687 GeV)2 and for "physical" m = 6.6 MeV with Dq = 128.0/(0.700 GeV)^. For 
all three curves, the chemical potential takes the same value /xcep = 164 MeV. 

We find that when the quark mass value varies from m — 5.2 MeV to m = 8.0 MeV, the 
most singular point of the quark number susceptibility emerges always along the /.fcEP = 
164 MeV line, while the corresponding Tcep has a very minor shift from 116 MeV to 118 MeV 
(the dashed and solid line, respectively). As we change the coupling parameter from its 
physical value in the present model Dq = 128.0/(0.687 GeV)^ to Dq = 128.0/(0.700 GeV)^ 
(with quark mass fixed at its physical value), the corresponding CEP moves slightly to 
TcEP — 114 MeV, with its /xcep remaining at 164 MeV. Lastly, as for the parameter A, 
the effect of its variation can be canceled by taking a larger numerical cutoff in practical 
calculations, since it appears in the exponential as the denominator. Also notable here is 
that the critical fine does not change much, either for m — 5.2 MeV, the critical chemical 
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potential at vanishing temperature is /Xc = 280 MeV, for m = 8.0 MeV, /Xc = 289 MeV, 
and for the new couphng parameter Do = 128.0/(0.700 GeV)^, fic = 281 MeV, while all 
the pseudo-critical temperatures take almost the same value as that for the "physical" 
parameters. 

The conclusion is that the location of the CEP (especially /icep) and the critical line is 
relatively robust with respect to the variations of the model parameters around their physical 
values. 



IV. SUMMARY AND CONCLUSIONS 

The primary goal of this paper is to identify and locate the CEP through a continuum 
study of the quark number susceptibility. To this end, we first derive a model-independent 
integral formula expressing the quark number susceptibility in terms of dressed quark propa- 
gator and dressed vector vertex. The latter objects are then confronted in the rainbow-ladder 
DSE-BSE framework with a confining separable model gluon propagator. With these results, 
the quark number susceptibility is calculated. The quark number susceptibility measures 
the response of the quark number density to changes in the baryon chemical potential and 
is of particular interest for exploring the CEP in that it is expected to diverge at the CEP. 
Our calculation reproduces this diverging phenomenon, which may be detectable in heavy 



ion collision experiments 



3, 



25| 



The location of the CEP our calculation predicts lies at a reasonable value. Interestingly 
the /icEP is coincident with recent experimental estimate and more acceptable than NJL-type 
model predictions. In this regard, we emphasize that the confining separable DSE model 
admits a momentum-dependent effective interaction, which is not the case for NJL-type 
models. Moreover, the rainbow-ladder truncation scheme for DSE-BSE respects the axial- 
vector Ward-Takahashi identity (AV-WTI) [sO^, which relates closely to the chiral symmetry 
and its dynamical breaking, and the usual Abelian WTI, which plays an essential part in 
our analysis in this paper. 



Nevertheless, our model study says nothing about the recently prevalent termino^ 



ogy m 



QCD thermodynamics, namely the strongly coupled Quark-Gluon Plasma(sQGP) 371]. rj/s, 
rather than the quark number susceptibility, is believed to provide a good characteristic of 
sQGP liquid and its surprisingly small value makes a current focus jisl . [s^], with which we 
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are also concerned. 
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